Instabilities in the dissolution of a porous matrix 
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. . . A reactive fluid dissolving the surrounding rock matrix can trigger an instability in the dissolution 

, front, leading to spontaneous formation of pronounced channels or wormholes. Theoretical inves- 

■ ligations of this instability have typically focused on a steadily propagating dissolution front that 
' separates regions of high and low porosity. In this paper we show that this is not the only possible 

dissolutional instability in porous rocks; there is another instability that operates instantaneously 
^—^ . on any initial porosity field, including an entirely uniform one. The relative importance of the two 

mechanisms depends on the ratio of the porosity increase to the initial porosity. We show that 
' the "inlet" instability is likely to be important in limestone formations where the initial porosity 

is small and there is the possibility of a large increase in permeability. In quartz-rich sandstones, 
0^ ' where the proportion of easily soluble material (e.g. carbonate cements) is small, the instability in 

' the steady-state equations is dominant. 

S: 

O ! I- INTRODUCTION 

CD . 

Dissolution is of fundamental importance in a variety of geological systems, including diagenesis, karst formation, 
Q ' aquifer evolution and melt migration 0]. It also plays an important role in a number of engineering applications, 
•i-H , such as dam stability fsl, CO2 sequestration risk assessment of contaminant migration in groundwater and 
stimulation of petroleum reservoirs [6] . The dynamics of a dissolution front in a porous matrix is complex, even under 
(-H , laminar flow conditions, with a number of possible feedback loops between reaction, diffusion and flow Q. These 
O i' feedback processes may trigger an instability in the front, leading to the formation of high permeability zones [sl-flO}. 
\ As a result of this reactive-inflitration instability, long, finger-like channels or "wormholes" are formed, which can 
^ 1: carry active reactant deep into the matrix [Tl| . 
^ Previous investigations of the reactive-infiltration instability (8l-[To| assumed that a steadily propagating dissolution 

front, separating regions of high and low porosity, forms first; subsequently an instability in this initially planar front 
T-H develops. However, we have found that dissolution may be unstable from the very beginning, even before a reaction 
00 , front is formed. Starting with an entirely uniform porous material, we have analyzed the competition between the 

■ development of a planar dissolution front and the growth of instabilities in that front. Here we focus on two limiting 
[ cases: first, the instability in the dissolution of an entirely homogeneous porous matrix (the "inlet" instability) and 
• second, the instability in a steadily propagating reaction front (the "front" instability). In the convection-dominated 
' limit the inlet instability shows a strong wavelength selection, similar to that found in dissolving fractures [l^ , whereas 
, the front instability is unstable over a broad range of wavelengths. However, the addition of a small diffusive flux 

^ ' changes the character of the front instability, introducing a strong wavelength selection in that case as well; this has 
been overlooked in previous work • The relative importance of the two mechanisms for unstable growth of the 

dissolution front depends on the amount of soluble material in the rock. 



X 



II. DISSOLUTION IN POROUS MEDIA 

We begin with standard equations for the dissolution of a porous matrix. The superficial fluid velocity in the pore 
space, V, is given by Darcy's law: 

v^-Ki<t>)^, (1) 
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FIG. 1: Sketch illustrating the two regimes of dissolution: the porosity field (j){x) during the initial dissolution t < tmax, Eq. ((9|, 
is show as the solid line and the steadily advancing porosity profile t ^ tmax (fT3)l. is shown as the dashed line. 



where the permeability K{(j)) is taken as a function of the fluid fraction (j) and /x is the fluid viscosity. The continuity 
equation takes the form 

dt^ + V-v^Q, (2) 

assuming the fluid filling the pore space is incompressible. 

Transport of reactants and products are described by the convection-diffusion equation; for simplicity we consider 
a single species with a dilute concentration field c, 

dt {(t>c) +v-Vc = V ■ D{(P, v)(f> ■ Vc - R{c), (3) 

where D{(j), v) accounts for dispersion of reactants by flow through the porous matrix, and R{c) describes the rate of 
consumption of reactants by the porous matrix. Finally, we have an equation for the dissolution of the porous matrix, 

VsoiCsoidt<i) = vR{c), (4) 

where Csoi is the concentration of the solid species, and Uaoi and v are stoichiometry numbers. 

The boundary conditions are constructed from the following considerations. We assume that the material is initially 
homogeneous with a porosity 0o and semi-infinite in extent, < a; < oo. Far from the inlet the material is undissolved 
so that Vx{x — >■ oo) = wqi where vq is a constant; in effect we fix the far-field pressure gradient to give the desired flow 
velocity. A constant concentration Cin is maintained at the inlet (x = 0). 

We will make a number of additional simplifications which do not affect the overall conclusions, but which allow 
for a simpler analysis and clearer presentation. First, we assume a linear kinetic equation, 

R{c) = -fcsc, (5) 

where k is the reaction rate, s is the specific surface area and c is the reactant concentration. In this equation we 
assume that the reaction rate is sufficiently small that diffusion within the pore spaces can be neglected. Otherwise a 
more complicated constitutive equation is required to take account of mass transport within the pores. The specific 
surface area can be related to the grain size by assuming the grains are spherical, in which case s = 6(1 — 0)/c? where 
d is the grain diameter. Although s increases as the grains dissolve, to keep the derivation simple we take s to be 
constant; then R{c) = —ks^c with sq the specific surface area in the initial state. We use the Carman-Kozeny equation, 
K = 0'^/5sg, to describe the permeability, maintaining the assumption that the specific surface area remains constant 
during dissolution. Finally, we will only include diffusion perpendicular to the flow, with a constant diffusivity D. It 
is possible to relax any or all of these assumptions at the cost of a more involved analysis. 



III. ONE-DIMENSIONAL ANALYSIS 



We first consider a one-dimensional solution of the dissolution equations ([I])-©, which will be used to establish a 
reference porosity profile for the linear stability analysis. Taking the flow to be in the a:-direction, a one-dimensional 
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porosity field evolves according to the following equations: 

dtcj) + dxVx = 0, 

dt{(t)c) + VxdxC = -fcsoc, (6) 

VsolCsoldt(f> = vksQC. 

The form of Eq. (jH) suggests that there are characteristic length and time scales in the dissolution of the porous 
matrix; the penetration length of the reactant concentration Ip = {kso/vo)~^ and the time to double the inlet porosity 
^2 = (t^sq/'/'o)"^! where 7 — vcin/ (vsoiCsoi)- From now on we scale length by Ip and time by t2- In addition we scale 
concentration by the inlet concentration, c — c/cim the superficial fluid velocity hy vq, i) = v /vo, and the porosity by 
00, 4> = 4>/4'o- Then Eq. ([6]) can be written in dimensionless form: 

-fdtcj) + dxVx = 0, 

-idt{(j>c) +yAc = -c, (7) 
dt(j) = c. 

The time scale for the evolution of the porosity is controlled by the coefficient 7, which is typically less than 10^^ 
in calcite dissolution, even for strong acids. Here we take the limit 7 — > 0, when the velocity and concentration fields 
become slaved to the porosity. The superficial fluid velocity is then constant throughout the domain, Vx = 1, and the 
concentration field is time independent, 



The porosity grows linearly in time. 



c{x) = e-^. (8) 



= 1 + ie-^, (9) 



until only insoluble material remains; i.e. when (f){x,t) = 4>max, the maximum porosity (relative to its initial value). 
Equation [9] remains valid for (dimensionless) times up to 

tmax — 4^max 1 — (10) 

at which point the reaction at the inlet stops. For times larger than tmax a propagating front develops, which moves 
through the porous matrix. Although the full solution is complicated, even in one-dimension, we can consider the 
long-time limit where a front propagates steadily with a constant velocity Vf. In this case we create a new coordinate 
system relative to the front position Xf — Vft, which is defined as the rightmost point where — 4>max, 

x'{x,t) — X — Vft. (11) 

The dimensionless dissolution equations in the new coordinate system are (still in the limit 7 — >■ 0): 

dx'Vx' = 0, 

Vx'dx'C = -c, (12) 
dt^ - Vfdx'cj) = c. 

At steady-state {dtcf) — 0) the porosity profile is again exponential, 

^{x') ^l + vf-^e-""' . (13) 
The boundary condition at the front <j){x' = 0) = 4>max determines the velocity of the front, 

«/ = A-i, (14) 

A steadily advancing porosity profile, Eq. (IT5|) . is sketched in Fig. [1] alongside the initial [t < tmax) profile, 
Eq. ([9]). Surprisingly, both profiles are unstable with respect to infinitesimal perturbations, although the nature of 
the instabilities is qualitatively different in the two cases. 

In our analysis we neglect any coupling between the development of the dissolution front and mechanical compaction. 
In a related problem of the channeling instability in an upwelling mantle ^ , mechanical compaction was shown to 
stabilize long-wavelength perturbations without affecting the value of the most unstable wavelength. 
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FIG. 2: Growth rates of the inlet instability. The solid line is for the reference state frozen at to = without diffusion [H = 0). 
The dashed lines show results for increasing H: top to bottom H — 0.1, H — 1, and H = 10. The dotted lines show results for 
the reference state frozen at later times: top to bottom to = 1 and to = 10 with H — 0. 



IV. LINEAR STABILITY ANALYSIS 




We develop a linear stability analysis by considering a small sinusoidal perturbation about the one-dimensional 
solution; 

4>{x, y, t) = ^r{x, t) + g{x) sin(uy)e'^* (15) 

where 0r is the one-dimensional reference porosity field given by Eq. ([9|) or Eq. (fT3| and g{x) is the amplitude of 
the perturbation; the dimensionless wavevector u = 27r/A, where A is the wavelength (in units of Ip). Note that the 
base state of the inlet instability, Eq. is itself time-dependent; here we use an approximation [l3[ in which the 
base state is frozen at a specific time, Iq, and the growth rate is then determined as if the base state were time- 
independent (the quasi-steady-state approximation). A solution of the inlet instability problem for the analogous case 
of a narrow fracture can be found in Szymczak and Ladd [l^ . In a similar fashion, after some straightforward but 
lengthy manipulations, we obtain a third order differential equation for the perturbation g, 

(16) 

<Pr / ' <Pr 

where (j)'^ = dx(t>r- The dimensionless constant 

H = Dkso<t)o/vl (17) 

can be written as the ratio of the Damkohlcr number Da — k/vo and the Peclet number Pe — vo/Dso<j)o, H = DaPe~^. 
Boundary conditions on g follow from the boundary conditions on the porosity, pressure, and concentration fields: 

g{x - 0) = 0, 

[5^e-5(a:)]^^„=0, (18) 
[d.e-g{x)]^^^ = 0, 

There are three independent solutions for g{x) but one of these is eliminated by the far- field boundary condition. Since 
the amplitude of g is arbitrary only one more condition is needed to fix the solution. The final boundary condition is 
then used to obtain the dispersion relation for the growth rate (t(m). Results for the growth rate of the inlet instability 
are shown in Fig. [5] for different values of H and to ■ 

The reference porosity field for a steadily- moving front, Eq. (jl3l) . is time independent in the coordinate system 
relative to the front. The resulting equation for the perturbation in porosity is fourth order. 




+ Hu'<Pr]e\a-Vfdx)g^^g. (19) 



In this case there are four independent solutions for g{x)^ but two of these are always eliminated by the far-field 
boundary condition. Care must be exercised in applying the remaining two boundary conditions at the front itself, 
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FIG. 3: Growth rates of the front instabihty. The sohd hne shows the case A —J- 0, without diffusion {H = 0). The dashed 
Unes indicate increasing diffusion: top to bottom H — O.OOf , and H — O.OI. The dotted lines show the effect of increasing A: 
top to bottom A = 2 and A = fO with H = 0. 



which is no longer planar but has small sinusoidal perturbations as well. As in the inlet instability, one boundary 
condition is sufficient to fix the solution to within a multiplicative constant, and the final boundary condition then 
provides the dispersion relation illustrated in Fig. [S] 



V. DISCUSSION 



The solid lines in Figs. [2] and [3] indicate the most rapid growth of each instability: the inlet instability at the onset 
of dissolution (ig 0) and the steady-state instability in the limit that the front develops instantaneously (A — >• 0). 
In both cases we are considering the convective limit H — 0. The inlet instability shows a strong wavelength selection, 
with a maximum growth rate at a wavenumber Uma-x ~ 1-3. Even when the flow rate is small {H > I) the growth rate 
is similar to the convection-dominated case, although the peak is pushed to longer wavelengths. On the other hand 
the front instability grows across a broad spectrum of wavelengths, limited only by the eventual onset of diffusional 
stabilization (dashed lines). Diffusion has a much larger effect on the front instability than on the inlet instability. 
Even for weak diffusion, H — 0.001, the shape of the dispersion curve changes qualitatively; short wavelengths 
are stabilized and a maximum growth rate appears. This suggests that the difFusionless limit considered in earlier 
work [To| is singular, and any amount of diffusion will lead to the appearance of a maximum in the dispersion 
curves. 

The solid curve in Fig. [2] was obtained by freezing the base state at to — 0. If the base state is frozen at a later time, 
then the growth rate is reduced as shown by the dotted curves. However, the peak growth rate remains at almost 
the same wavelength, independent of io, so a single mode {u ~ Umax) dominates until the onset of non-linear growth. 
In the case of the front instability, the growth rate is sensitive to the porosity contrast across the dissolution front 
(Fig.[3|). For small A, the growth rate increases monotonically towards a limiting value, in agreement with Hinch and 
Bhatt [9]. However, as A increases, the dispersion curve flattens, and at around A « 3.2 a maximum appears even 
in the convective limit, indicating a similar competition between reaction and convection as in the inlet instability. 
Interestingly, the most unstable wavelength in the front instability approaches that for the inlet instability when 
A > 1. 

The analysis outlined in Sec. II VI includes two parameters, the ratio of Damkohler and Peclet numbers, H = DaPe~^, 
and the porosity contrast, A = (t>max — 1- Fracture dissolution, where an initially small (less than 1mm) aperture 
opens essentially without limit, can be considered as the limiting case A oo. In this case a steadily advancing 
front never develops and Eq. (fTB)) applies. We can obtain the corresponding equations for fracture dissolution by 
mapping the porosity (j) onto the dimensionless fracture aperture h = h/ho; the initial fracture aperture ho is related 
to the specific surface area by ho = 2/so- For a typical calcite fracture - reaction rate fc = 2.5 x 10^^ cms^^, aperture 
ho = 0.2 mm and hydraulic gradient of 1% [l3| - the dissolution is convection dominated {H ~ 10"^) and the growth 
rate of the instability follows the solid line in Fig. [2] Thus, for these parameters we would expect to see fingers with 
an initial spacing of the order of lirlp / Umax ~ 1 m, but the competition for flow between growing wormholes will 
eliminate many small channels and the spacing observed at later times will be larger. We have previously suggested 



that this instability in the fracture dissolution front initiates the formation of underground caves in limestone [12 1. 

A convection-dominated instability can also develop in highly porous and permeable limestone formations. Crinoidal 
limestone, for example, can have a speciflc surface area as low as 100 cm~^ and a permeability as high as 10~* cm^ [l5| . 
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For a gravitationally-driven flow, H ^ 1Q~^ and the penetration length is of the order of 1 cm. The corresponding 
wavelengths of an unstable front are of the order of several centimeters for the inlet instability and one or two 
centimeters for the front instability. The inlet instability is active from the beginning of the dissolution process with 
fingers developing on a dimensionless time scale t ^ 1. At the same time there is uniform dissolution as indicated 
by Eq. (|9]), which in the absence of the inlet instability would establish a dissolution front on a dimensionless time 
scale i A, Eq. (fTO|) . The relative importance of the two instabilities will therefore depend on the porosity contrast 
A, which typically ranges from 1 (0o = 0.5) to 10 (00 — 0.09). For large values of A, we would expect substantial 
fingering from the inlet, but in more porous formations a steadily advancing front would appear first. 

At the opposite extreme from fracture dissolution is the limit A — > where the porosity contrast is very small, 
although there can still be a substantial change in permeability. This limit includes sandstones, where the inert quartz 
grains are interspersed with a soluble cement [TgI . I17| . In such cases a steadily advancing front will always form prior 
to the instability. 

For smaller grain sizes (sq ~ lO^'cm"^) the flow rate is extremely small {H ^ 1) and a depletion layer develops 
behind a steadily propagating front Q . The concentration field in the region behind the front is significantly perturbed 
by diffusion of solute from the reaction zone and we we can no longer assume a constant concentration at the front. 
In this diffusive limit a uniform front is stable at low flow rates ^8|], but above a critical flow rate, fingers are predicted 
to develop with a wavelength of order D/vq. 

Experimental observations of dissolution patterns in porous media @, 0, [HI are typically made under conditions 
corresponding to reservoir acidization, so the reaction rates are much higher than is typical under geophysical condi- 
tions. For example, when acidizing calcite with hydrochloric acid the reaction rate is of the order of O.lcms"^ @. 
Taking sq ~ lO'^cm"^, corresponding to micron-sized grains, and 7 = 0.01, the front is established on a time scale 
of (ksQ-f)^^ w 0.1s. Similar considerations apply to other experiments and indeed to acidization in the field. In 
acidization it seems clear that the instability develops out of a steadily propagating front [1, [l^ . We reiterate that 
diffusion cannot be ignored in the acidization instability, even under high flow rate conditions (Fig. [3]). 
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